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Abstract —The task of approximation of points with circu¬ 
lar arcs is performed in many applications, such as polyline 
compression, noise filtering, and feature recognition. However, 
development of algorithms that perform a significant amount of 
circular arcs fitting require an efficient way of fitting circular arcs 
with complexity 0(1). The elegant solution to this task based on 
an eigenvector problem for a square nonsymmetrical matrix is 
described in 0. For the compression algorithm described in Q, 
it is necessary to solve this task when two points on the arc are 
known. This paper describes a different approach to efficiently 
fitting the arcs and solves the task when one or two points are 
known. 

Index Terms —arc fitting; optimization; compression; general¬ 
ization 

I. Introduction 

The purpose of this paper is to solve the task of fitting a 
circular arc to a set of points (or segments) with complexity 
0(1) when two points on the arc and moments up to the fourth 
order are known. 

In papers Q and Q, fitting a circle is done by finding such 
circle with center {xc,yc) and radius r, which minimizes the 
next equation 

n 2 

+ {yi-ycf^ , ( 1 ) 

i=l 

where {xi,yi) are Ath point, i = l..n. 

This formula is minimizing square differences between 
square distances from the circle center to the points and square 
of the radius. The solution is found by using only the moments 
of {xi,yi) with complexity 0(1). However, this leads to bias 
in the estimation of parameters |[^ see pp. 368-370]. Suppose 
that each point has been fitted with error. Substituting it in 
(fT) gives 

n 2 

X • 

After simplifying this equation 

n n 

X (2^ ■ -f . 
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Assuming that are small compared to radius r and 
neglecting higher orders 

n 

^ 2 ) 

i=l 


From this formula, it is clear that the fitting is trying to 
decrease the radius to minimize the equation 0. When the 
points cover only a small part of a circle, the estimated center 
of the circle can move toward the arc to reduce the radius. 
Despite that, the errors are increased, while overall penalty 
0 is decreased. The smaller the angle of the arc, the worse 
the effect. Fig. [T] shows the case when fitting by moments 
produces a circle with too small a radius. 



Fig. 1. Comparison of fitting an arc using different approaches. The 
comparison is performed for the arc with 72°, and uniform noise in the circle 
is proportional to 10 percent of the ai‘c radius. A total of 1, 000 random points 
were simulated along the arc with uniform steps. The black arc is a ground 
truth arc. The black dots are source points. The red circle is a solution based 
on fitting squares of distances (see ^). The green circle is solution based 
on fitting distances (see j^). The blue circle is a solution described in this 
paper (appr oximate solution of |5| found by one iteration of the algorithm 
described in|Appendix B: Minimization of Multidimensional Function f (a:),| 
\x £ - 

Minimization of the next equation was suggested in 0: 















It has this error formula; 


This covers all possible values of {xc,yc,T) and gives a 
significant advantage for finding the minimum by removing 
(4) the third and fourth order variables in the numerator of (|^: 
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E 


i=l 


This doesn’t have a problem like in (|^. However, mini¬ 
mizing (|^ requires an iterative approach, which analyzes all 
points {xi,yi) leading to an algorithm with complexity 0(n). 
The efficient algorithm to find the minimum of ([^ is described 
in ||6). 


II. Algorithm 

The solution to remove the square root from was 
developed in 0’ ® see p. 675], 0, and 0- From 0. 
there comes an idea that dividing Q by 4r^ and minimizing 
it will produce a result closer to 0 because it is close to 
0. The approximation based on the Taylor expansion of the 
square root by the first two terms gives exactly this solution. 
Approximation of y/x at a; = 1: 

y/x I + - {x - 1) + O = 2 2 ' 

Applying this approximation to 0 



Unlike 0, this formula can be minimized using only 
moments. The direct solution, based on conformal geometric 
algebra, is described in ||T]. The solution is based on finding 
eigenvalues of a square nonsymmetric matrix |[T] see (24)]. 
This can be done by Schur factorization. The eigenvector 
corresponding to the minimal non-negative eigenvalue is the 
solution. Care needs to be taken for the cases when there is a 
solution with close to zero eigenvalue due to round-off error. 

In this paper, another approach to minimization of Q will 
be considered. 


III. Minimization of 0 

Suppose we have a good estimate {xe,ye,i’e)- Let’s opti¬ 
mize it in the next form: 

(^Xe + Ax, ye + Ay, -f Ax^ + Ay^ -f Ar^ . 


/(Ax, Ay, Ar) = 
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+ + yl - rl) + 

-|-2Ax (Xe — Xi) + 
+2Ay (ye - yi) - 
—ArJ 


4 (r^ + Ax^ + Ay2 -|- Ar) 


Writing it from moments 


/(Ax, Ay, Ar) = 

V + Vx ■ Ax + Vy • Ay + Vr • Ar-f 
+'Wa:,£c • Ax^ + Vy^y ■ Ay^ + Ar^-f 
+Vx,y ■ Ax ■ Ay -f Vx,r • Ax • Ar + Vy^r ■ Ay ■ Ar 
4 (r| -I- Ax^ -I- Ay2 + Ar) 


where 

V = (M4_o + 2Af2,2 + -^ 0 , 4 ) ~ 

— 4 (M3_o + Mi^ 2) Xe — 4: (M2,1 + Mq^s) ye + 

+ SMi^l • Xe • ye + 2 AT 2 P • Zx + 2Mq^2 ‘ Zy — 

— 4 (Ml 0 • Xe + Mq i ■ ye) z + , 

Vx = 4:(— (Ma^o + Mi_2) + (3M2_o + Mo,2) Xe + 

+2Mi^l ■ ye - 2Mo,l • Xe • ye - Mifi ■ Zx+Xe- z) , 

Vy = 4{— (M2_i + Mq^s) + (M 2 P + 3Mo^ 2) J/e + 

-f2Mi_l • Xe - 2Mi_o ■ Xe-Ve- Mq,! ■ Zy + ye ■ z) , 

Vr = —2 (M 2 P + Mg,2 ~ 2 {Mifi ■ Xe + Mg,! • ye) + z) , 
Vx,x = 4 {M2,0 - 2Mi_g • Xe + Xg) , 

Vy,y = 4 (Mg _2 — 2Mo^i • ye -\- ye) i 

Vx,y = 8 (Mi_i - Mg^i • Xe - Mi_o ■ ye + Xe ■ Ve) , 

Vx,r = 4(Mi^g - Xe) , 

Vy^r = 4 (Mg,i - ye) , 

2,2 2 
z = Xe +ye -r^, 

Zx = 3Xe + 2/e ~ 
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Zy Xg -|- 3yg Tg, 
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Minimization of 0 can be done using the approach de- 
scribed in [Appendix B: Minimization of Multidimensional 
Function f (x), x € To use that approach, it is sufficient 






















to know the matrix of second derivatives up to the constant 
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where 


dx,x — 2^ {v Vx,x ‘ ‘^e) * 


dx, y — ^x,y ‘ 

dy, y = -{v - Vy,y ■ rl) -rl, 

dx, T — ( '^x “t” '^x,r ' ^e) * 

dy, r = {~'^y + • Tg) • Tg, 

dr,r =2 [v - Vr ■ rl + r®) . 

and the equation for directional search by direction 

i^rxxi rxy j Oi'Y'^ is 


V [Vx ■ ax + Vy ■ ay + Vr ■ ar) t+ 

{^XjX ‘ a^ “t“ '^y,y ' ^y ‘ ax ‘ Cty-t- 

~\~Vx * ax * 0^7' ~1~ Vy^r * ay * f 

4 (r^ + i + (ctx + cty) 

Looking at the numerator of it would be reasonable to 
take a solution of ([T]), described in |^, as a starting point. 

Only a few iterations are needed to converge beyond ma¬ 
chine precision. However, this approach is only approximate, 
and there is no need for such precision. In practice, one 
iteration is sufficient to get a good approximation. 

When estimation of the center is known, the best estimation 
of the radius can be easily found. However, it does not give 
any improvement in speed. 

An approximation of the sum of square deviations from the 
polyline to an arc with the center {xe,ye) and radius rg (see 
([^) is found from ^ by setting Ax, Ay, and Ar to zero and 
multiplying by n. 


4r2 


(7) 


The algorithm described in Q has an advantage for finding 
the global optimum, while the algorithm described in this 
paper can go to the local optimum. This is likely to happen 
when the arc is close to the line. Otherwise, the results are 
identical. 

The advantage of the approach described in this paper is the 
ability to reduce the amount of calculation by approximating 
the solution. 

I have implemented both approaches. Intel Math Kernel 
Library 11.2 was used to solve the nonsymmetric eigenvector 
problem in |jT). The approach described in this paper is several 
times faster. However, it is difficult to make a fair comparison 
due to the different ways of implementing and optimizing the 


code. When speed is not a concern, the approach described in 
Q is preferred. 

IV. Optimal Arc when One Point is Known 

For some tasks, one point on the arc is known in ad¬ 
vance. The arc should pass through that point. Knowing 
th e position of the center determines the radius: r = 

J (^{Xr - Xaf + iyc - yaf^, where {xa,ya) is a point on 
the arc. Substituting this into Q 

- Xcf + iyi - ycf^ - 
(^{Xc - Xaf + {yc - yaf'^ 



4 (^{Xc - Xaf + iVc - Vaf 


( 8 ) 


The solution described in section [11^ [Minimization of 
can be applied. The differences are that the optimization is 
performed in two-dimensional space and the initial solution 
can be found from the least squares approach. 

Multiplying the numerator and denominator of (|8]l by 
and replacing Xc ■ s and y^- shy Vx and Vy, respectively, and 


setting V = 


where 
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B= 0 


Xn 


i -Va 

-Xa -ya Xl + ylj 


© is a generalized Rayleigh quotient. Note that A and 
B are symmetric non-negative matrices. The solution can be 
found by the generalized singular value decomposition of ^fA 
and V®B. Square root matrices can be found from singular 
value decompositions such as 


= LlXj, 

Vb = LlXl, 


( 10 ) 

















where La and Lb are eigenvalue matrices of A and B, 
respectively; Xa and are eigenvector matrices of A and 
B, respectively. 

From generalized singular value decomposition for a/A and 
^fB follows 

^/A = UDa{Q,R)Q^. 

Vb = VDb {0,R) q ^, 

where U, V, and Q are orthogonal matrices; R is an upper 
triangular matrix. 

From ( [Tol l and o follows 

A = Q{0,RfD\ (0,i?) Q^, 

B = Q{0,RfDl (0,i?) Q^. 


The smallest ratio of squares of eigenvalues is the solution. 

/Vx Vy\ 

The center of the arc is recovered from u as ( —, — ■ 

Vs s / 


s s 

V. Optimal Arc when Two Points are Known 


There are tasks when the starting and ending points of 
the arc are known (or any two points lying on the arc). 
Therefore, the center of the arc should lie on some line of 
{xc, Vc) = (xp + a ■ t,yp + 13 ■ t), where {xp, yp) is a point on 
the line, {a, (3) is the direction of the line {a^ + /3^ = 1), and 
t is any value. K nowing the position of the center determines 

the radius: r = {{xp + a ■ t) - Xaf + {{yp A l3-t)- yaf, 
where {xa,ya) is one of the points on the arc. Substituting 
this in (j^ 


n ( {{Xi - {Xp + a ■ t)f + {y^ - {yp + P ■ t)f^ 

\ - {{{Xp + a ■ t) - Xaf + {{yp + (3 ■ t) - yaf"j ^ 

4 f ((xp + a-t)- Xaf + {{yp + (3-t)- yaf) 

( 12 ) 

It is easy to see that the minimum of can be easily 
found because it has the form of ([T3]l in [Appendix A: Finding 


the Global Minimum of the Ratio of Quadratic Equations An 
iterative algorithm is not needed to solve this problem. 

The algorithm described in Q finds an optimal polyline 
within the tolerance of the source polyline, with the minimum 
number of vertices, and among them, with the minimum sum 
of the square deviations from the optimal polyline. Extending 
this algorithm to support arcs requires efficient fitting of the 
arc from the known start and end points and evaluation of the 
sum of the square deviations from the source polyline to an 
arc. An approximate solution (0 can be used instead of direct 
evaluation 


VI. Example: Recovering Arcs in a Cadastral 
Dataset 

The approach described in this paper for efficiently fitting 
circular arcs is used in a compression algorithm, when ver¬ 
tices of the source polylines are not allowed to move. The 
algorithm minimizes the weighted number of segments (with 
penalty 2) and arcs (with penalty 3) while satisfying tolerance 
restrictions. Among all possible solutions, the solution with 


the minimum sum of square deviations is chosen. A dynamic 
programming approach was used to find the optimal solution 

ID- 

Algorithms with a very similar framework were described 

in 1^, (ni. 

Eitting of arcs in pO) was performed by checking tolerance 
when starting and ending vertices are fixed. It has the advan¬ 
tage of always finding an arc within tolerance; however, the 
computational complexity for each fitting is 0{n). This paper 
uses approximation to least squares fitting with complexity 
(3(1) described in section |v| [Optimal Arc when Two Points 
[are Knowml Although checking for the tolerance and proper 
sequence (zigzag) |T^ has complexity 0{n), it is only 

performed for optimal fits. 

The penalty function in o is a combination of perceptual 
and fitting errors. The perceptual error is 

d - sin-, 

where S is the segmentation penalty and a is the angle between 
adjacent segments. This gives preference to solutions with 
acute angles. 

An example is shown in Eig. |D The original arcs were lost 
due to digitalization, limitations of the format, projection, and 
so forth. The restoration of arcs is an important task because 
it is the original representation. Restoring original arcs creates 
cleaner databases and simplifies future editing. 



Fig. 2. Part of a parcel map with lost circular arcs. A compression algorithm 
was applied to this data. The Black lines are the source polylines, the red 
circles are vertices of the source polylines, and the green asterisks are resultant 
vertices. All original arcs were reconstructed. 


The optimizations used for fitting segments are described 
in ID- Initial fitting of circular arcs to different parts of the 
polyline was performed to reduce searching for arcs in a dy¬ 
namic programming approach. The processing speed is about 
0.03 seconds per 100 points on an Intel Xeon CPU E5-2670. 















VII. Future Work 
Described at the end of section |VT] [Example: Recovering! 


Arcs in a Cadastral Dataset reducing the search for arcs 


is only approximate and in some cases might skip possible 
solutions. There is an exact solution to check if an arc exists 
that covers all vertices within the specihed tolerance described 
in m section 7.4]. The solution is found using the closest 
and the farthest point Voronoi diagrams. The center of an arc 
corresponding to the minimum width covering all vertices is 
either a vertex of the closest or the farthest point Voronoi 
diagram or a point on the edge of the closest and the farthest 
point Voronoi diagrams Gil p. 167]. The closest and the 
farthest point Voronoi diagrams are dual with the closest 
and the farthest point Delaunay triangulations, respectively 
G 3 - Note that the farthest point Voronoi diagram includes 
only vertices on the convex hull. The algorithm to construct 
the farthest point Delaunay triangulation is similar to the 
algorithm for the closest point Delaunay triangulation and can 
be found in Another algorithm to construct the closest 
and the farthest point Delaunay triangulations for vertices on 
the convex hull is by mapping each vertex {x, y) to a vertex 
in three dimensional space {x,y,x^ + 2 /^) and constmcting a 
convex hull for them G 3 - 


VIII. Conclusion 

This paper describes an efficient approximation for fitting 
circular arcs. While all formulas are for a two-dimensional 
case, the algorithm can be generalized for higher dimensions 
(for example, htting a sphere to points). 

The direct solution to ht arcs is described in ||T]. This paper 
extends the solution to cases when one or two points on the 
arc are known. 

Because the solution is based on fourth orders, it has a 
negative impact on the precision of calculations. This can be 
solved by shifting data to the origin of a coordinate system 
and/or using floating point numbers with a larger mantissa. 

There is no evaluation of how well the ht is done. An 
additional algorithm is necessary to perform this check, as 
described in in see section 3]. 


The hrst derivative of ( [T3| l equals 

(ai • 6o - oo • 6i) + 2 (02 ■ bo - ao ■ 62) • x+ 

+ (02 • 61 - oi • 62)■ 

-9-■ ( 16 ) 

{bo + bi ■ X + b2 ■ x"^) 

From ( fTSj l, it follows that the denominator ([T^ is always 
positive in Q. Therefore, it is sufficient to work with the 
numerator: 

Co -f Cl • X + C2 • X^, ( 17 ) 

where cq = oi • hg — oq • bi, ci = 2 (02 ■ bo — ao ■ 62 ), 
and C 2 = 02 • — fli • 62 . 

From ( [l4| ), it follows that ( [T3| ) is not negative in Q. If the 
denominator of ( [T3] l has real roots, then ( [T3] l, when x is 
approaching any root, goes to +00 in Q and —00 in the 
complement of Q excluding roots (see example in Fig.l^l. 
Local extrema are found from roots of 03 (Figs'- There 
is a special case, when in ( [T3] l the numerator is equal to 
zero at one of the roots of the denominator. In this case, 
( pj] ) simplihes to the ratio of linear equations and doesn’t 
have any global minimum. 



Fig. 3. Example of G3- The area outside domain Q is shown in gray. 
Local extrema are shown by red circles, found as the solution of (Tt) (see 
Fig.|4). 


Appendix A: Finding the Global Minimum of the 
Ratio of Quadratic Equations 

ap + ai ■ X + a2 ■ x^ 

bp + bi ■ X + b2 ■ x^ ’ 

where Ui and bi are known coefficients i = 0..2. 

Coefficients should satisfy 

Vx, ao -I- Cl • X + 02 • x^ > 0. (14) 

Two cases will be analyzed separately: 

7. 62 7 ^ 0. 

The domain will be restricted to 

Q = {bo + bi ■ X + b2 ■ x^ > 0 } . (15) 

Notice, limits ( |T3] l when x —7 —00 and x —7 +00 are the 
same. 



Fig. 4. Example of ([T 3 corresponding to the function shown in Fig. 
Roots are shown by red circles. 

The global minimum can be found from roots of the 
quadratic equation ( [T7] i: 

a. If C 2 > 0 and the largest root of ( [T7| ) belongs to Q, 
then it is a global minimum. 

b. If C 2 < 0 and the smallest root of ([T^ belongs to Q, 
then it is a global minimum. 
















c. If C 2 = 0, Cl > 0 and the single root of ( |T7] i belongs 
to Q, then it is a global minimum. 

d. Otherwise, no global minimum exists. 

To summarize, the global minimum of ( [T3] l can be found 
by the next equation if the value inside Q 

if C2 = 0 A Cl >0, 

if C 2 7 ^ 0 A D > 0 A Cl <0, 

— — if C 2 7 ^ 0 A D > 0 A Cl = 0, 

C2 

if C 2 7 ^ 0 A D > 0 A Cl >0, 
otherwise, 

(18) 

where D = c\ — 4co • C 2 is discriminant of ([T7J. 

2. 62 = 0. It is sufficient to evaluate the solution of the next 
equation to show that this case can be properly solved by 

m 

Oq + Oi • a; + 02 • 

X 

where Oi are known coefficients i = 0 .. 2 . 

The domain will be restricted to 


' Co 
Cl 

\/D - Cl 

2 c 2 

sign (C 2 ) • ^ 
2 Co 

y/~D + Cl 

no solution 


Q = {a: > 0} . 

The first derivative multiplied by equals 

—oo + 02 • a;^. 


From that global minimum 



no solution 


if Oo • 02 > 0 , 
otherwise, 


(19) 


Notice that solution ([T^ is equal to solution ( [T8] l. There¬ 
fore, it is sufficient to use for both cases. 

Another way to prove that the solution for the case 
gives the proper solution (when 62 = 0 ) is to consider 
lim of (fTSll. 

62^0 ’—' 


Appendix B; Minimization of Multidimensional 
Function f (a;), x e K” 

Suppose that the minimum of f {x) along any direction can 
be easily found as well as second derivatives. 

The next algorithm is suggested; 

a. Let i = 0. Define the starting point Xq. 

b. Find the second derivative matrix at Xi and find all 
eigenvectors. 

c. For each eigenvector, from Xi point, search along the 
eigenvector direction for minimum Xi+i. Set i = i + 1. 
Because the number of eigenvectors is n, this step in¬ 
creases the index of x by n. 

d. If Xi is close to the minimum with enough precision (for 
example, by comparing with the previous estimate Xi_„), 
then stop; otherwise, go to step 


In the case of quadratic functions, this algorithm converges 

to the minimum in one iteration consisting of searching from 

any starting point by n direction. 
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